Redox Entropy of Plastocyanin: Developing a Microscopic View of 

Mesoscopic Polar Solvation 
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We report applications of analytical formalisms and Molecular Dynamics (MD) simulations to 
the calculation of redox entropy of plastocyanin metalloprotein in aqueous solution. The goal of 
our analysis is to establish critical components of the theory required to describe polar solvation 
at the mesoscopic scale. The analytical techniques include a microscopic formalism based on struc- 
ture factors of the solvent dipolar orientations and density and continuum dielectric theories. The 
microscopic theory employs the atomistic structure of the protein with force-field atomic charges 
and solvent structure factors obtained from separate MD simulations of the homogeneous solvent. 
The MD simulations provide linear response solvation free energies and reorganization energies of 
electron transfer in the temperature range 280-310 K. We found that continuum models universally 
underestimate solvation entropies, and a more favorable agreement is reported between the micro- 
scopic calculations and MD simulations. The analysis of simulations also suggests that difficulties 
of extending standard formalisms to protein solvation are related to the inhomogeneous structure 
of the solvation shell at the protein-water interface combining islands of highly structured water 
around ionized residues along with partial dewetting of hydrophobic patches. Quantitative theo- 
ries of electrostatic protein hydration need to incorporate realistic density profile of water at the 
protein-water interface. 



I. INTRODUCTION 

Calculations of the thermodynamics of hydrated 
biopolymers present a significant challenge to theoretical 
algorithms. In many cases the problem can be treated by 
numerical simulations with force fields assigned to the 
biomolccule (solute) and water (solvent) fi The obvious 
difficulty is the large computational load and still existing 
uncertainties in the treatment of the long-range electro- 
statics. The problem, however, becomes more nontrivial 
when derivatives of thermodynamic potentials, e. g. re- 
dox entropy, need to be computed or when the solvation 
thermodynamics changes on the time- and length-scale 
unattainable to standard Molecular Dynamics (MD) pro- 
tocols, c. g. in problems related to protein foldingi^i^ In 
all such cases, coarse graining of the system is required 
and that can be done on various length-scales 4*^ Dielec- 
tric continuum algorithms, solving the boundary Pois- 
son problem, are computationally very efficient. In this 
approximation, all length-scales below the largest dis- 
tance of microscopic correlations (density and/or polar- 
ization) arc averaged out into a continuum surrounding 
the solute. These approaches are normally represented 
by either direct solution of the Poisson-Boltzmann equa- 
tion on the real-space grid^ or even more approximate 
formalisms under the umbrella of the generalized Born 
approximatiouj^ 

When the cavity cut out by the solute from the con- 
tinuum dielectric is properly parametrized, equations of 
continuum electrostatics provide a reasonable estimate of 
the solvation Gibbs cnerg y^'^i-'^" The fundamental prob- 
lem of this approach is that the local structure of the 
solvent around the solute, averaged out in the contin- 
uum limit, is what effectively forms the dielectric cavity. 
While this structure can be parametrized by choosing 



proper van der Waals (vdW) radii,— this parametriza- 
tion needs to be re-done every time the thermodynamic 
state of the solvent changes. This difficulty makes contin- 
uum formalisms unreliable for the calculation of deriva- 
tives of the Gibbs energy, for instance the entropy of 
solvatioufi^iHiiiii^ In addition, the surface of a protein is 
highly heterogeneous combining hydrophobic patches ex- 
posing non-polar residues and hydrophilic patches made 
of ionized/polar residues. While the water structure is 
rigid around ionized residues, probably resembling the 
well-studied case of solvation around simple ions, water 
is much less structured at hydrophobic patches with the 
potential for dewettingi^ or/and oscillations of the wa- 
ter occupation! ^^1^^ It is clear that simplistic continuum 
does not represent this complex reality;^ and one needs 
to incorporate the ability of the solvent to fluctuate into 
the solvation model. 

The goal of this paper is to extend the microscopic view 
of solvation in polar solvents, which we have been devel- 
oping in the past in application to small and medium-size 
solutes r^ii^ii^ to solvation of solutes of mesoscopic di- 
mension, biopolymers in the first place. The length-scale 
of this problem presents an obvious obstacle to numeri- 
cal simulation techniques. On the other hand, the same 
length-scale allows one to hope that some of the short- 
range features of the solvent structure around the solute, 
making solvation of small molecules so specific, might av- 
erage out on a larger scale. If this averaging is realized for 
solvation of biopolymers, it would allow coarse-grained 
models to efficiently operate in this field complementing 
direct numerical simulations. Our approach to the prob- 
lem is to coarse-grain the solvent response into a number 
of solvent correlation functions (structure factors) rep- 
resenting the nuclear modes affecting electrostatic sol- 
vation. The microscopic nature of the solvent response 



is then incorporated into the wave-vector dependence of 
these structure factors efficiently filtering out the length- 
scales insignificant for solvation. 

This study poses the central question for the future 
development of such techniques: What are the solvent 
modes which play the central role in the thermodynam- 
ics of mesoscopic polar solvation and what are the the- 
ory ingredients critical for capturing the basic physics of 
large-scale solvation? We study this problem here by 
carrying out extensive MD simulations of solvation of 
plastocyanin (PC) in TIP3P water in the temperature 
range 280-310 K. This fully atomistic approach is com- 
pared to continuum electrostatics and to our microscopic 
algorithm, operating with k-space correlation functions, 
which was designed to scale efficiently on the mesoscopic 
length-scale. 

Plastocyanin from spinach is a single polypeptide chain 
of 99 residues forming a /3-sandwich, with a single cop- 
per ion coordinated by 2 sulfurs from cysteine and me- 
thionine and 2 nitrogens from histidine residues (Fig. [1]). 
The presence of the copper ion, which can change redox 
state, allows PC to function as a mobile electron carrier 
in the photosynthetic apparatus of plants and bacteria. 
It accepts an electron from ferrocytochrome / and diffu- 
sionally carries it to another docking location where the 
electron is donated to the oxidized form of Photosystem 
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FIG. 2: Distribution of the positive and negative charge on 
the surface of the protein. The positively and negatively 
charged residues are shown, respectively, in red and blue. The 
copper ion is shown in green. 



lysine residues with amino groups protonated (Fig. ^. 
The asymmetric charge distribution located on the pro- 
tein surface creates the dipole moment of 2200 D in the 
oxidized state (Ox) and of 2470 D in the reduced state 
(Red) , both numbers are calculated relative to the center 
of partial charges. 

The redox potential of the protein includes a compo- 
nent from the local ligand field of the active site and 
the Gibbs energy of solvation. The computation of the 
former requires quantum mechanics, making the prob- 
lem of calculating the overall redox potential a very non- 
trivial exercise i^^'^°i'^^ Calculations of solvation thermo- 
dynamics can be reasonably accomplished using partial 
atomic charges parametrized from quantum calculations 
in the vacuum. The experimental input comes from mea- 
surements of redox entropy^ii^i^^ since the temperature- 
independent ligand-field component is expected to vanish 
in the temperature derivative. 



II. MICROSCOPIC SOLVATION MODEL 



FIG. 1: Structure of plastocyanin: the active site includes 
copper ion (green), 2 histidines (blue), methionine (red), and 
cysteine (orange) residues. 

The redox thermodynamics of PC has been char- 
acterized expcrimentallj«2ii22ii^ and combined quan- 
tum/simulation calculations have been done as well.— i^ 
The early focus of the theoretical studies had been on un- 
usually high redox potentials of copper proteins, which 
was assigned to the non-traditional distorted tetrahedral 
coordination on the copper ion.— i^ii^ In particular, the 
Cu-S bond to methionine is unusually long and is actu- 
ally broken in the reduced state of PC at pH < 3. 8^2^ 
The protein is also highly charged at pll~ 7 (—9.0 in re- 
duced state and —8.0 in oxidized states). The charge is 
made by 15 negatively charged deprotonated residues (9 
glutamic and 6 aspartic acids) and six positively charged 



The principal idea of the microscopic solvation model 
is to reduce the problem of solvation of an arbitrary so- 
lute in a polar solvent to a formalism combining two 
major blocks: electrostatics of an isolated solute and 
non-local correlation functions of the pure solvent. The 
idea of assembling separate solute and solvent proper- 
ties in a solvation model is obviously not new going 
back to Born^ and Onsager— and all the subsequent 
development of continuum electrostatics in application 
to solvatioui^iiiSiiSi^ The advantage of our approach is 
in avoiding the necessity to know the microscopic solute- 
solvent structure, which is the main complexity of mi- 
croscopic solvation models and is also their main advan- 
tage when the problem is successfully resolved by ei- 
ther solving integral equations^ or by applying time- 
dependent^i^ or equilibrium^i^i^ density functional 
methods. Inserting a solute into a dense liquid creates a 



significant distortion of its structure, and the incomplete 
account of the couphng between the short-range density 
profile around the solute with the long-range polarization 
field is perhaps the weakest part of our formulation when 
applied to small solutes4^ [This deficiency is almost com- 
pletely off-set by averaging of the density profile around 
a nano-scale solute (see below).] On the other hand, a 
strong side of our formalism, its ability to treat solvation 
of large solutes of irregular shape and arbitrary charge 
distribution , ^' becomes particularly useful in applica- 
tion to protein solvation. 

The reduction of the many-body solvation problem to 
an irreducible representation in terms of a few basic cor- 
relation functions depends on the symmetry of the solute- 
solvent interaction potential. The number of correlation 
functions is known to grow with increasing the rank of the 
solvent multipole included in the interaction potential.'^^ 
Solvent dipoles are for the most part sufficient for sol- 
vation in polar liquids (^ in which case the solute-solvent 
interaction potential Vqs ( "0" and "s" are used for the so- 
lute and solvent, respectively) is a sum of pairwise inter- 
actions of the solute electric field Eo(r) with the solvent 
dipoles 



Vas = - 



N 

E 



m. 



Eo(r,). 



(1) 



Here, m'- is the dipole moment characteristic of the bulk 
state of the solvent; m' is usually higher that the vacuum 
dipole m because of the collective field of the induced sol- 
vent dipoles4^ For instance, the dipole moment of water 
in the liquid state, 2.4-2.6 D^ is higher than the gas- 
phase dipole of 1.83 D. 

We will focus on the electrostatic component of the 
chemical potential of solvation fiQ^ which contains all the 
information relevant to electrostatic solvation. Linear re- 
sponse approximation (LRA) significantly simplifies the 
problem and provides several equivalent routes to fiQs- 
One can consider the full interaction between the atomic 
charges of the solute and the solvent and determine /ios 
as half of the average solute-solvent interaction energy^ 
fJ-Os = (Vos)/2. Alternatively, one can use the second 
cumulants, {{SVos)^) or (((5Vbs)^)o'^ In the first cumu- 
lant, the angular brackets (...) refer to an ensemble av- 
erage over the fluctuations SVqs in the solvent in equilib- 
rium with the full charge distribution of the solute. For 
the second cumulant, (. . ■ )o implies that all the charges 
of the solute have been set to zero, and fluctuations of 
the solvent in the solute vicinity are regulated only by 
short-range solute-solvent interactions, molecular repul- 
sions in the first place. In the LRA, the two cumulants 
are equal,— which physically means that the solute elec- 
trostatic forces do not significantly change the solvent 
structure around the solute established by the prevalence 
of short-range repulsionsi^ Computer simulations for the 
most part support this picture^i^i^ with a few excep- 
tions of very strong solute-solvent electrostatic coupling 
found for small solutesJ^i^i^ 



This observation opens up a significant simplification 
of the calculation algorithms. Instead of solving the in- 
homogeneous problem of restructuring the solvent in an 
external field of the solute, it appears to be sufficient to 
look at the statistics of solvent fluctuations around the 
repulsive core of the solute. This strategy is used here 
and we will base our calculations on the relation 



/io,s = i/3/2){{SVQsf)o, 



(2) 



where 6Vos = Vos - {Vos)o and (3 - 1/ikBT). 

By using the interaction potential according to Eq. ((T|), 
one can re-write Eq. ^ in the form typical for Gaussian 
(LRA) models of solvation^i^ 

- /ios = ^Eo(ki) * x(ki,k2) * Eo(k2). (3) 

Here, the 2-rank tensor x0^ii^2) is the response 
function^ of the system composed of a dipolar solvent 
and a solute to a weak field of the solute. The inhomo- 
geneous character of the problem is reflected by the fact 
that x(ki,k2) depends on two wave-vectors, ki and k2, 
separately and not on ki — k2 , as is the case with response 
functions of homogeneous solvents. The asterisk in Eq. 
Q refers to both the tensor contraction and integration 
in inverted k-space. In addition, Eo(k) is the Fourier 
transform of the electric fleld of the solute deflned by the 
integral limited to the solvent volume fl: 



Eo(k) 



Eo(r)e*''-''dr 



(4) 



The shape of the solute thus enters both the response 
function x(ki, k2) and the fleld Fourier transform Eo(k). 
The charge distribution of the solute, which determines 
the electric field Eo(k), is given by its electronic density 
and is commonly represented by partial atomic charges. 
The main challenge of this formalism, as well as of 
other Gaussian solvation theories)^ is how to connect 
the inhomogeneous response function ;;\;(ki,k2) to the 
shape of the solute repulsive core and the self-correlation 
functions of the solvent modes affecting solvation. Two 
modes naturally appear in most theories: dipolar (orien- 
tational) polarization and density fluctuationsj ^^i^°'^^i^° 
For the former, the combination of axial symmetry in- 
troduced by the wave-vector k with the vector character 
of the dipolar polarization P(k) allows one to split the 
2-rank tensor Xs(k) = (/3/rj)(|(5P(k)p) into the longitu- 
dinal and transverse dyads i^ 



x.(k) = I [J^^^fc) 



J^S^ik)] , 



(5) 



where 3^ — kk, 3^ = 1 — kk. In Eq. ([5]), y is the effective 
density of both permanent and induced dipoles in the liq- 
uid which commonly appears in theories of dielectrics*^ 



In Eq. g 
particle. 



y = (47r/9)/3?7i'V + (47r/3)pa. (6) 

a is the dipolar polarizability of the solvent 



The scalar functions S^{k) and S'^{k) in Eq. ^ are, 
correspondingly, the longitudinal and transverse struc- 
ture factors of dipolar fluctuations of the homogeneous 
solvent (see below) . The k = values of these structure 
factors are related to the dielectric constant, e^, by the 
following equations: 



5^(0) = (6, - l)/{3yes), 
S^iO) = (6, - l)/(3y). 



(7) 



Also, the trace of Xs(0) over the Cartesian projections is 
proportional to the Kirkwood g-factor— 



5K = J[^^(0) + 25^(0)]. 



(8) 



The expansion of the solvation chemical potential in 
the Mayer functions corresponding to the solutc-solvcnt 
interaction potential leads to the following form for the 
response functionHiiS 



where 



X(ki,k2)=xp(ki,k2)+Xd(ki,k2), (9) 



Xp(ki,k2) = x,s(ki)(5ki,k2- (10) 



In Eq. (Uni), f^(5ki,k2 = (27r)-^(5(ki - k2) is the Kronecker 
symbol and Xd(ki,k2) in Eq. ([9]) is the component of 
the response originating from the local fluctuations of 
the solvent density around the solutoi^ 



Xd(ki,k2) = (32//8^)(l - S{ki))0o (ki 



(11) 



Here, S{k) ~ N^^{\Sp{'k)\'^) is the density-density struc- 
ture factor of the homogeneous solvent and N is the num- 
ber of solvent molecules. In addition, 0o(k) is the Fourier 
transform of the step function 6*0 (r) defining the solute 
shape. It is equal to unity for r inside the solute and is 
zero otherwise. 

The problem with the direct perturbation result in Eq. 
([Q]) is that it contains the transverse polarization response 
function oc S'^{k) diverging in its continuum, fc — > 0, limit 
as the solvent dielectric constant goes to infinity [Eq. ([7])] • 
The problem is really caused by the non-spherical shape 
of the solute. The electric field of the solute charges is 
longitudinal. However, when the symmetry of the solute 
is different from the symmetry of the charge distribution 
in a sense that the cavity boundary does not coincide 
with the equipotential surface, the Fourier integral in Eq. 
dl]) generates a transverse component of Eo(k). Notice 
that this is always the case when electron transfer reac- 
tions are considered!^ A transverse component in Eo(k) 
then results in a "transverse catastrophe" for solvents of 
high polarity. The problem was well recognized in early 
studiesi^i^i^ which suggested to use only the longitudi- 
nal component of the field Eo(k). As a matter of fact, the 
problem lies in the response function of the dipolar po- 
larization field which needs to be re-normalized with the 
account of the solute repulsive core, a procedure similar 



to applying boundary conditions to the Poisson equation 
of continuum electrostatics. 

The Li-Kardar-Chandle r^^i^'^ Gaussian model allows 
one to achieve a correct renormalization of the inhomo- 
gencous polarization response function Xp(ki,k2) elimi- 
nating the "transverse catastrophe" i^ This approach in- 
troduces another simplification by replacing all the short- 
range solute-solvent interactions by hard-core repulsions. 
This simplification, however, leads to an exact solution 
for the k-space response functions with the result)^ 

Xp(ki,k2) = Xs(ki)<5k„k. -X"(ki)0o(ki-k2)xs(k2). 

(12) 
The second summand in Eq. (fT2|l is the correction of the 
response function of the homogeneous solvent, appearing 
in Eq. ((T0|) . by the repulsive core of the solute. The 
response function x"(Mi) then incorporates both Xs and 
the information about the solute shapej ^'^'^^ 

A direct substitution of Eq. P^ into Eq. ^ results 
in a 6D integral convolution in k-space which is not nu- 
merically tractable. In order to arrive at a computation- 
ally efficient procedure, a mean-field approximation was 
introduced)^ which replaces the inhomogeneous electric 
field of the solvent inside the solute by a mean cavity 
field: 



F„^X 



Stt 



Eq • Dr^T, 



where 



/ 



2(e. - 1) 



2e, 



1 



(13) 



(14) 



Here, Dr ~ 3ff — 1 is the 2-rank dipolar tensor with 
f = r/r. Fq becomes the Onsager reaction field^ for a 
spherical solute with point dipole located at the center. 
The mean-field approximation reduces the problem of 
calculating the solvation thermodynamics to a numeri- 
cally tractable 3D integral in k-space. The chemical po- 
tential of solvation then becomes a sum of two compo- 
nents arising from the longitudinal (L) and transverse 
(T) polarization fluctuations, fi^^ , and a third compo- 
nent arising from the density fluctuations, /Iq^: 



MOs 



Mo. 



T 



Mo. 



(15) 



The transverse component fiQ^ is defined by k-integral 
of the transverse projection of the solute filed, Sj (k), 
with the transverse polarization structure factor 



^^L^9KS'^m 



dk 



,3y 

8n J (27r)3 



lEiiWS'ik). (16) 



The transverse field component is defined by subtracting 
the longitudinal projection 



(17) 



E^(k) =k(k-E 

from the total inverted-space field Eg: 

Ej(k)=Eo(k)-E^(k). 



(18) 



Equation (J16p is the main result of the apphcation 
of the Gaussian model^ to polar solvation. It replaces 
S'^{k) of the direct perturbation expansion in Eqs. ([9|) 
and (1121) with the renormalized function: 



3yS^{k) ^ i3yS'^{0)/gK)S^ik). 



(19) 



The fc = limit of the transverse response changes from 
3yS'^{0) = e^ - 1 to 3(es - l)/(2e, + 1) thus eliminat- 
ing the "transverse catastrophe" of direct perturbation 
expansions. 

The component fi^^ of the solvation chemical potential 
in Eq. (fT5| is obtained by inverted-space integration with 
the longitudinal polarization structure factor: 



L 



dk 



3y 

Sir J (27r)' 



5^fc) 



T|2 _ 



E^\'~\Ei\'f 



^ 



^0 



(20) 

There is a significant physics behind the appearance of 
the transverse field in the brackets of Eq. ([^0]) . Longitu- 
dinal dipolar polarization is short-ranged and thus does 
not propagate over macroscopic distances. On the con- 
trary, transverse polarization is long-ranged. Therefore, 
inducing transverse polarization modifies the electric field 
acting on the solvent dipoles resulting in the second term 
in the brackets in Eq. pO)) . 
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FIG. 3: Diagram of the computational algorithm. 

The density component in Eq. (jlSp can formally be 
obtained by multiplying the response function XdO^i, ^2) 
[Eq. pT|) ] with the Fourier transforms of the electric field 
and integrating over ki and k2. This, however, results 
in a 6D convolution integral to be avoided in numerical 
applications. An alternative approach is to use direct- 
space integration when the density component becomes 



^^os 



= 3y^(r)*J-i[(l-5(fc))0o(k)] 



Here, F{r) ~ (87r)^^i?Q(r) is the density of the electro- 
static field energy and the asterisk indicates integration 
in real space over the volume Q occupied by the solvent. 
In addition, !F~^ is the inverse Fourier transform of the 
function indicated in the brackets. 

Solvation by the overall dipolar polarization, includ- 
ing nuclear and electronic components, was considered 
in the formalism outlined above. For solvation prob- 
lems relevant to spectroscopy and charge-transfer reac- 
tions nuclear component of polarization needs to be ex- 
tracted. This is achieved by replacing the density y [Eq. 
(I6|)] of all, permanent and induced, dipoles in the equa- 
tions above with the density of permanent dipoles only, 
y ^ Up = {4:n/9)Pm'^p. In addition, the k = values 
of the structure factors need to be modified to account 
for screening of the dipolar interactions by the high- 
frequency dielectric constant eao- The fc = values for 
these nuclear structure factors, S^''^{k), now become^^ 



5n(0) 






^e;')/{3y,), 
eoo)/(32/p). 



(22) 



Once the A; = values for the structure factors are 
fixed by Eqs. ^ and ^^, the scalar functions S^{k) 
and S {k) can be calculated from our parametriza- 
tion scheme, parametrized polarization structure factors 
(PPSF).— This analytical route to the polarization struc- 
ture factors is tested here by comparing the results of 
solvation calculations employing the PPSF to the direct 
use of S'"'^{k) from MD simulations (see below). 



III. COMPUTATIONAL ALGORITHM 



The computational algorithm is outlined in Fig. [31 The 
solute is parametrized by coordinates r^, vdW radii Qj, 
and partial charges qj of the atoms. The electric field of 
the solute is calculated at points r„ of the N x N x N 
grid built on the L x L x L cube: 



N„ 



Eo(r„)=^ 



*■(!■« -r^) 



(23) 



i=i 



(21) 



where Nq is the number of solute charges. The ar- 
ray Eo(r„) is converted to inverted space by using fast 
Fourier transform techniquci^ The field Eo(k) is split 
into longitudinal and transverse components and used in 
k-integration in Eqs. pB|) and ([^0)1 with the correspond- 
ing structure factors of the dipolar polarization. As is 
illustrated in Fig. \3\ the calculation input is subdivided 
into two separate components related to the solute and 
solvent properties. Details of calculation for each of these 
are given below. 



A. Solute 

The Fourier transform of the Coulomb field is condi- 
tionally convergent. Therefore, in order to avoid numcr- 



ical divergence in the Fourier integral, real space is di- 
vided into three regions: hard core of the solute (region 
1), region outside a sphere of radius R (region 3), and the 
region between the solute surface and the sphere (region 
2, Fig. 2]). The Fourier integral then becomes 



(24) 



Eo(k)= / Eo(r)e-*'^-'-dr + E;,, 

Jr<R 

where the first integral is taken over region 2 and 



Efl(k) = 



Eo(r)e 



— zk-j 



dr. 



(25) 



r>B. 



The center of the sphere is taken at the center of the 
charge distribution defined by the relation 






(26) 



The radius of the cut-off sphere R is chosen to minimize 
the part of the grid which is used in numerical calcula- 
tion of the Fourier transform. In our calculations, the 
radius R is chosen by adding the solvent diameter a to 
the largest distance from r^ to the solvent-accessible sur- 
face (SAS) of the solute (vdW radii of the surface atoms 
plus the radius of the solvent molecule) . 




FIG. 4: Separation of real space into regions for the calcula- 
tion of the Fourier transform of the solute electric field [Eq. 
(|24[l ]. The Fourier transform is calculated numerically in re- 
gion 2 and analytically [Eq. (|27|l ] in region 3. The field is 
set equal to zero within the hard repulsive core of the solute 
(region 1). 

The Fourier transform outside the sphere can be eval- 
uated analytically. For the location of charges relative to 



the center of charge given as s^ 



To, the solution 



for Efl can be obtained by expanding Eo(r) in Sj/R < 1: 

-ijkR) 

^ (27) 



E^(k)=-47re^g,£(^ 



Jn- 



s,p;:_i(cos0,)-kp;,(cos0,) 



Here, cos Oj = Sj -k, j„ {x) is the spherical Bessel function, 
and Pn{cos9j) is the Legendre polynomial. 



B. Charging Scheme 



The formal charge of the copper ion is +2 and of the 
cysteine sulfur is —1 in the oxidized state of PC. The 
charge is, however, delocalized among the ligands and 
the metal center. The main factor in this delocaliza- 
tion is a strong covalency of a copper-sulfur (Cys) bond. 
Calculations by Solomon and co-workers^ assign 40% 
of spin density of an unpaired electron to copper and 
36% to cystein's sulfur in the ground state of oxidized 
protein. The extent of delocalization varies significantly 
depending on the level of quantum mechanical calcula- 
tions uscd]^i^i2£ The electron-nuclear double resonance 
(ENDOR) experiments^^ which require additional cali- 
bration on quantum calculations, result in the following 
net charges on the residues coordinating copper>22 —0.25 
(His), -0.51 (Cys), and -0.04 (Met). The more recent 
mapping of the electron spin density to NMR relaxation^ 
gives a much lower extent of delocalization: —0.11 (Cys), 
-0.025 (His), (Met). 

The uncertainties in the extent of electron delocaliza- 
tion pose the question of their impact on the calcula- 
tion of the redox thermodynamics. In order to study this 
question, we have performed calculations of the solvation 
part of the redox potential and the corresponding entropy 
using different charge sets. Set I is chemically fake assum- 
ing charge +2 on copper in Ox state and the net charge 
of —1 on cysteine. The negative charge is placed on cys- 
teines sulfur in addition to —0.23 from CHARMm22 pro- 
tein paramctrization. The rest of the protein charges are 
from the standard CHARMm paramctrization. The re- 
duced state for Set I is obtained by changing the metal 
charge to -|-1. The charges for copper and its four lig- 
ands are summarized in Table [H Set II is based on the 
charging scheme listed by UUmann et al^ for the oxi- 
dized state of PC. The reduced form is obtained by plac- 
ing an extra negative charge on copper and its three lig- 
ands, Ns (His87), Ns (His37), S'^(Cys84), and 5.^(Met92) 
in proportion extracted from NMR experiments (Table 
IJ),^ Finally, a third charge distribution is completely 
parametrized at the DFT level for the charges and force 
constants of the copper and ligand atoms and consistent 
with the Amber force fieldi^ In addition. Amber FF03 
parametrization2i was applied to all non-ligand residues 
(Set II). There were various numbers of TIP3P water 
molecules for each of the charge distributions: 5,874 (Set 
I), 5,886 (Set II), and 4,628 (Set HI). 

We ran separate simulations (ca. 5 ns) for each charg- 
ing scheme to find that the results are not strongly af- 
fected by the choice of atomic charges (Table |TT| . This 
was also noticed in some other recent simulationSi ^^i^^ 
We have therefore implemented charge scheme II in all 
simulations reported here since it presents a reasonable 
balance between being simple and realistic. 



TABLE I: Atomic partial charges for copper and its four ligands in the reduced (Red) and oxidized (Ox) states of PC. 



Set 


Cu 


Nr 


Red 


S,^ 


c d 
^5 


Cu 


Nr 


Ox 


S^^ 


c d 

as 


I 
II 


1.0 
-0.49 


-0.7 
-0.445 


-0.7 
-0.495 


-1.23 
-0.369 


-0.09 
-0.24 


2.0 
0.35 


-0.7 
-0.42 


-0.7 
-0.47 


-1.23 
-0.26 


-0.09 
-0.24 


"HisS? 
''His37 

'^Met92 























TABLE II: Temperature dependent (Vba)ox (eV) for PC(Ox). The results are obtained from MD simulations, DelPhi calcula- 
tions (with vdW and SS cavities) , and from NRFT calculations. The calculations were done with three charge distributions of 
the active site (I-III, see text for description). 







MD 






DelPhi(vdW) 






DelPhi(SS) 




NRFT" 


T/K 


I 


II 


III 


I 


II 


III 


I 


II 


III 


II 


280 


-68.29 


-69.46 


-68.37 


-104.86 


-102.07 


-99.76 


-46.90 


-47.47 


-46.37 


-104.6(-37.2) 


285 




-69.69 




-104.81 


-102.03 


-99.71 


-46.88 


-47.45 


-46.36 


-103.3(-36.5) 


290 


-70.73 


-66.01 




-104.76 


-101.98 


-99.67 


-46.87 


-47.43 


-46.34 


-102.1(-35.8) 


295 




-66.97 




-104.71 


-101.93 


-99.62 


-46.84 


-47.41 


-46.32 


-100.7(-35.1) 


300 


-67.06 


-66.09 


-65.84 


-104.65 


-101.88 


-99.57 


-46.81 


-47.39 


-46.30 


-99.3(-34.5) 


305 




-68.58 




-104.59 


-101.82 


-99.52 


-46.80 


-47.37 


-46.28 


-98.3(-33.9) 


310 




-67.51 


-66.68 


-104.53 


-101.77 


-99.46 


-46.79 


-47.35 


-46.26 


-97.1(-33.3) 



"The numbor8 in the parenthese8 indicate the density component 
of the equilibrium interaction energy, (Vbs)Q^. 



C. Solvent 

The polarization structure factors entering the equa- 
tions for the solvation chemical potential are character- 
istics of the homogeneous solvent. They can be obtained 
numerically by averaging the projections of dipole mo- 
ments Gj on an arbitrary chosen direction of the k- vector, 
k = k/fc: 






iVi-Tij 



(28) 



where r^ = r^ — Vj and N is the number of liquid dipoles. 
Unfortunately, experiment does not provide spatially 
resolved correlators of dipoles in polar liquids and one has 
to resort to using either computer simulations or liquid- 
state theories. Parametrizing homogeneous structure fac- 
tors by computer simulations is a very attractive avenue 
for studying hydration because of continuously improv- 
ing empirical potentials for water— on one hand and the 
reliance of biological applications on aqueous solvation 
on the other. For the problems related to derivatives 
of thermodynamic potentials, e. g. entropy and volume 
of solvation, the structure factors need to be tabulated 
at different temperatures and/or pressures. In this pa- 
per, polarization structure factors S'^'^(fc) were obtained 




k/A 



FIG. 5: Longitudinal (L) and transverse (T) polarization 
structure factors of TIP3P water [Eq. (I28p ] calculated at dif- 
ferent temperatures from ~ 25 ns MD trajectories. Also 
shown is the PPSF calculation at T = 300 K. 



from ~ 25 ns MD trajectories of TIP3P water— at dif- 
ferent temperatures (Fig.E]). 

In parallel to simulations, we have used dipolar struc- 
ture factors from our PPSF parametrization schemeJ^ 
This approach is based on the analytical solution of the 
mean-spherical approximation for the fluid of dipolar 
hard spheres^^ which is parametrized to give fc = values 
from Eqs. ^ and (g^jjliii The PPSF parametrization 
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gives solvation free energies essentially identical to those 
obtained with the structure factors from simulations, as 
was also found for a smaller polypeptide solute in our 
previous publication.— It appears therefore that the lo- 
cal tetrahedral order of water, which is of course not cap- 
tured by dipolar hard spheres, is not significant for the 
energetics of polar solvation dominated by orientational 
correlations ruled by dipole-dipolc forces. The density 
fluctuations are, on the other hand, dominated by repul- 
sions. We used, therefore, the density structure factor 
from the Percus-Yevick solution^ for a hard-sphere fluid 
of the same density as water in calculations of the density 
response in Eq. (PT|) . 

We need to emphasize here that using dielectric con- 
stants of model dipolar fluids would give us wrong results. 
The dielectric constant of a molecular liquid is affected 
by short range molecular correlations through the Kirk- 
wood factor, hydrogen bonds and molecular quadrupoles 
arc among significant factors4^ We account for all these 
effects in the PPSF scheme by using the experimental, 
either from computer or laboratory data, dielectric con- 
stants in Eqs. ([7]) and ^^. Once the fc = limit is set 
up by the experimental input, the behavior of S^''^{k) in 
the range of fc < 2TT/a {a is the solvent diameter) is well 
reproduced by solutions obtained for model dipolar flu- 
ids. These solutions will fail at fc > 27r/a, where a is the 
characteristic distance between partial charges within the 
solvent molecule. This part of the spectrum of polariza- 
tion fluctuations is correctly captured by models based 
on interaction-site integral equations^ but that range 
of wave-vectors normally does not contribute to the sol- 
vation energy. In fact, the range of fc-values relevant for 
the solvation problem is limited by fc < 2'k/R, where R is 
the characteristic dimension of the solute. For large so- 
lutes, only the long-wavelength part of the polarization 
structure factors is really needed for the solvation energy 
calculations. As is shown in Fig. [5l there is a mismatch 
between the PPSF longitudinal structure factor and MD 
simulations. However, this difference makes no effect on 
the calculated solvation energies. 

We can summarize our results on parametrizing the 
solvent properties by stating that the model fluid of dipo- 
lar hard spheres can serve as a reliable reference system 
for calculations of polar solvation given the macroscopic 
properties, the density of dipoles y and the dielectric 
constant Cs, have been taken from experiment (either 
laboratory or computer). The theory thus adds an ad- 
ditional parameter y to the dielectric constant used in 
electrostatic solvation theories to produce a fully micro- 
scopic solvent response. In practical applications of the 
theory (e.g. in case of solvation in ambient water pre- 
sented below), the parameter y needs to be calculated 
from the molecular properties of the solvent. We use the 
1-R Wertheim theory^S to calculate the effective dipole 
moment of the solvent (see Ref. |8l| for comparison to 
simulations). The solvent input is thus made by five 
parameters: {cr, p, m, a, Cg}. One needs in addition the 
high-frequency dielectric constant Coo for the reorganiza- 



tion energy calculations and the temperature slopes of 
two dielectric constants, as well as the isobaric expansiv- 
ity, for the solvation entropy calculations. The big ad- 
vantage of the PPSF scheme is that all these parameters 
have been tabulated for many solvents commonly used in 
solution chemistry making our method broadly applica- 
ble to solvation calculations in polar molecular solvents. 
Despite the fact that the dielectric constant is sensi- 
tive to local correlations, the polarization structure fac- 
tors in the long-wavelength limit are fully determined by 
dipolar correlations general for all polar liquids and not 
much sensitive to details of the local structure which is of 
course very different in water than in a hard-sphere dipo- 
lar fluid. There are several advantages to using dipolar 
hard spheres as the reference system. First, all thermo- 
dynamic and structural properties are controlled by only 
two parameters, the reduced density pa^ and the dipolar 
density y. Second, this system is well characterized both 
analytically and numerically. It has served many times as 
a starting point for developing theories of polar liquids^ 
similarly to the role played by the fluid of hard spheres 
in theories of non-polar liquidsi^ Once that stated, we 
however want to stress that the theory itself is based on 
the structure factors of an arbitrary polar medium with 
the Gaussian fluctuation spectrum and is not limited to 
a choice of any particular reference system. 



IV. SIMULATIONS PROTOCOL 

Amber 8.0^ was used for all MD simulations. The ini- 
tial configuration of PC was created using a protonated 
version of the X-ray crystal structure at 1.7 A resolution 
(PDB: lag62^). This initial configuration of the protein 
was first minimized in vacuum by the conjugate gradi- 
ent method for 10,000 steps to allow the protein to re- 
move any bad initial contacts. Then the system was sol- 
vated in a rectilinear box with several thousand TIP3P 
molecules )2^ providing at least two-three solvation shells 
around the protein. To neutralize the charge, a num- 
ber of sodium ions equal to the total charge of the pro- 
tein were added. The protein was then relaxed for a few 
thousand steps while water and sodium were position- 
ally constrained. Finally, the entire system containing 
solvent, counterions, and protein was energy minimized 
in 100,000 steps. 

Next, the system was heated in a NVT ensemble for 30 
ps from K to the desired temperature followed by vol- 
ume expansion in a 1 ns NPT run. NVT production runs, 
following density equilibration, lasted from 6 to 18 ns. 
The last 5-10 ns at the end of each trajectory were used 
to calculate the averages. The timestep for all MD simu- 
lations was 2 fs, and SHAKE was employed to constrain 
bonds to hydrogen atoms. Constant pressure and tem- 
perature simulations employed Berendsen barostat and 
thermostat, respectivelyi^ The long-range electrostatic 
interactions were handled using a smooth particle mesh 
Ewald summation with a 9 A limit in the direct space 



sum. The total charge for the protein was —9.0 for the 
reduced state and —8.0 for the oxidized state. 



V. RESULTS 

The calculations presented here are focused on two 
properties: the solvent portion of the redox chemical po- 
tential, A^s, and the solvent reorganization energy As, 
both corresponding to the half reaction 



PC(Ox)^ 



PC(Rcd)' 



(29) 



The former can in principle be calculated as the differ- 
ence of solvation chemical potentials in the Red and Ox 
states. However, this approach involves calculating the 
difference in two large numbers, which is computationally 
unreliable. Instead, wc use the linear response approxi- 
mation to calculate A/i^ according to the equation: 



Here, Eq = (E?'' 



Red 



M?/ = -AEo*x*Eo. 



(30) 



E;,'''=^)/2 and AEo 



ERcd 



e: 



Ox 



are 



the mean and the difference of the electric fields in the 
Red and Ox states. Similarly, the solvent reorganization 
energy is calculated from 



As = -AEo * X * AEq. 



(31) 



Equation pip applies to the reorganization energy of 
non-polarizable solvents employed in computer simula- 
tions. For laboratory data, nuclear polarization should 
be separated from the overall solvent polarization and the 
response function x is replaced by the nuclear response 
function Xn as explained above and in more detail in Ref. 



A. Redox Thermodynamics 

The solvation thermodynamics calculated here can be 
related to experimental redox entropies reported by mea- 
suring the temperature dependence of the standard or 
midpoint electrode potentials i^^'^^i'^^ An electrochemical 
experiment corresponds to bringing a solution containing 
given numbers of oxidized and reduced reagents, which 
are not necessarily in equilibrium (the ratio of their num- 
bers is not a Boltzmann factor), in contact with a metal 
electrode. The equilibrium is established between the 
electronic subsystem of the redox pair and the electrode 
in such a way that the electrode is charged and its elec- 
trochemical potential fi is shifted from the vacuum Fermi 
energy ep by the electrostatic potential 0: ^ = tp — ecf). 

The numbers of the oxidized and reduced forms of 
the redox pair, A^ox and A^Red, are assumed to be large 
enough so that they are not affected by charging the 
electrode. The electrochemical potential of the electrode 



than becomes equal to the absolute electrochemical po- 
tential of the redox couple in the solution^^^ The lat- 
ter can be found from simple statistical arguments. The 
grand-canonical free energy of two fermionic subsystems 
of A^ox and A^Rod electronic levels is^^ 



pn = -A^ox In (l + e'3('^-'ox)A 
-A^R,dln(l + e''('^-^--)), 



(32) 



where eox and eRod are the average energies of the elec- 
tronic levels in the corresponding redox states. The 
chemical potential is then found by requiring that the 
derivative —{d^/dii)T is equal to the total number of 
electrons A^Rod- For the energy gap between Ox and Red 
states greater than k^T , this requirement results in the 
Nernst equation^ 



M 



eOx + ERcd 



fcBrin(Arox/AfRcd), (33) 



in which the standard potential is given by the mean of 
the average electronic energies 



eox + fiRcd 
Ye ' 



(34) 



The same result follows from the use of the stationary 
condition (zero electrode current) for the rates of reduc- 
tion and oxidatior 



,88 



^OxCOx = fcRcdCRod, 



(35) 



where Cox/Rcd are the surface concentrations. By using 
the Marcus equation for the reaction rate^^ 



fcox/Rod ex exp 



-P 



(eOx/Rod - m) 



2l 



4As 



(36) 



and neglecting the logarithmic correction including the 
ratio of two surface concentrations, one gets equal rates 
at /x ~ (eox + eRcd)/2. The double- well Marcus free en- 
ergy surface for the electrode electron transfer is then 
symmetrical as illustrated in Fig. [6l This picture bears a 
clear similarity with the formation of the Fermi level in 
the forbidden band of a semiconductor, as was noticed 
by Reissi^ 

The electronic energies are given by the sums of their 
vacuum components, Cox/Rcd' ^'^'^ ^^*^ interaction of the 
electric field of the electron Eg with the polarization of 
the solvent in equilibrium with the total electric field of 
the molecule in the solutionSL 



eox/Rod - eox/Rcd - Ee * X * Eq 



Ox/Red 



Taking into account that 



E, 



ERcd 




EOx 




AEo, 



(37) 



(38) 



one gets from Eqs. ([34)) . ((37)) . and (|38p the commonly 
used connection between the standard electrode potential 
and the solvation part of the redox free energy 



^Rod 



■^Ox 



A^s 



2e 



(39) 
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TABLE III: Redox thermodynamics of PC (eV). 







A/i,(TIP3P) 






A,(TIP3P) 






A.(H20) 




T/K 


MD 


NRFT 


vdW" 


SAS" 


MD 


NRFT 


vdW^ 


SAS 


NRFT" 


vdW 


SAS 


280 


-3.05 


-6.21 


-9.63 


-7.15 


0.845 


1.37(0.59) 


3.645 


0.690 


0.82 


2.728 


0.417 


285 


-2.73 


-6.15 


-9.62 


-7.14 


0.814 


1.35(0.58) 


3.641 


0.685 


0.80 


2.725 


0.417 


290 


-2.71 


-6.10 


-9.62 


-7.14 


0.646 


1.33(0.57) 


3.636 


0.685 


0.79 


2.722 


0.417 


295 


-3.14 


-6.04 


-9.61 


-7.14 


0.565 


1.32(0.56) 


3.630 


0.685 


0.78 


2.718 


0.417 


300 


-2.82 


-6.01 


-9.60 


-7.14 


0.435 


1.30(0.55) 


3.623 


0.684 


0.77 


2.714 


0.417 


305 


-2.94 


-5.94 


-9.59 


-7.13 


0.426 


1.29(0.54) 


3.620 


0.684 


0.76 


2.710 


0.417 


310 


-2.53 


-5.89 


-9.59 


-7.13 


0.543 


1.27(0.53) 


3.618 


0.684 


0.74 


2.710 


0.416 



"Poisson-Boltzmann calculations with vdW radii assigned to the 
protein atoms (standard vdW cavity). 

''Poisson-Boltzmann calculation with solvent radius added to the 
radii of the protein atoms exposed to the solvent (solvent-accessible 
cavity, SAS). 

"^The continuum calculations are done with the dielectric constant 
of TIP3P water e^ = 95 and e^o = 1.0. The dielectric constant of 
the protein interior was put equal to unity in order to be consistent 
with the microscopic calculations. The temperature variation of the 
dielectric constant oide^/dT = —0.654 K~^ (Ref. |14|) was adopted 
for the entropy calculations listed in Tab. IIVI 

''Calculations in ambient water using too = 1.78, tj, = 78, 
deJdT = -0.398 K"!, and de^a/dT = -2.75 X 10"* K~^ . In 
addition, temperature expansion was included with the constant- 
pressure expansivity coefficient ap = 2.6 X lO"** K~^. 



electrode 



F(X) 




^l-<'=Red> 



FIG. 6; Contact of a redox pair with the metal electrode, 
eox — M fi'iid ^Rcd — M show the fluctuating energy gaps for 
reduction and oxidation electron transfer, respectively. The 
equilibrium electrochemical potential of the electrode is es- 
tablished when the equilibrium energy gaps are equal for the 
reduction and oxidation reactions [Eq. (|33p ]. The Marcus 
electron transfer parabolas, shown by the dependence of free 
energy G{X) on the energy gap coordinate X = eox — A*, are 
symmetric in this case producing equal oxidation and reduc- 
tion currents [Eq. (|35|l ]. 



tween experiments done in buffered protein solution a^^i^^ 
and our calculations at zero ionic strength. 

From Eqs. ([M|) . ([571) . and ([551) O'^ci can directly derive 
the equation for the solvation redox free energy 



A^s = ((AVbs)ox + (Ayos)Rcd) /2, 



(41) 



where AVqs is the difference in the solute-solvent interac- 
tion energies in the Red and Ox states and the averages 
are taken over the corresponding ensembles. The same 
average vertical gaps can be used to calculate the reor- 
ganization energy as 



A, = ((A-t^os)ox-(AK)s)Rcd)/2. 



(42) 



We need to caution here is that while Eq. ([341) is a 
statistical-mechanical result, Eqs. ([37l) -([42l) are based on 
the LRA for the solute-solvent interaction energy and 
might be affected by deviations from this approximation. 



B. Solute-solvent average energy 



where A^s is given by Eq. (|30p. The first term in this 
equation disappears in the temperature derivative re- 
ported experimentally2ii^i2^ 



dT 



As, 



„Rcd 



f^^^i. 



dT 



(40) 



We need to stress here that redox entropies in polar so- 
lutions are sensitive to the presence of electrolyte i^i2^ 
One therefore can expect only a qualitative agreement be- 



In addition to our NRFT formalism, we have used the 
dielectric continuum approximation implemented in the 
Delphi program suited in the solvation calculations. Di- 
electric constant of ambient water was used for the sol- 
vent continuum and Cg = 1 for the protein. This latter 
choice was driven by our desire to compare continuum 
and microscopic calculations of solvation thermodynam- 
ics since the latter does not assume any polarization of 
the protein. Table |TT] hsts the results of NRFT and Del- 
Phi calculations of the average energy (Vbs)ox of PC in 
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FIG. 7: Average solute-solvent interaction energy (Vbs)ox 
obtained from MD simulations (closed circles), NRFT (di- 
amonds), and DelPhi continuum calculations (vdW cavity, 
triangles). The closed diamonds refer to the total average 
energy including the polarization and density components, 
while the open diamonds denote the polarization component 
only. The dashed lines represent linear regressions through 
the points. 



4 

3 

"b 

1 





3 

"9 2 
1 




Ionized Surface Residues 
ASP9-0. 



81 
LYS77-N, 



;i, J i>i||i yi i ) ii. > ! im. ' .r.!,! ) i.iii. ') .. ' !. ' ..,: ! lmi..J.. 



Non-Polar Surface Residues 

PR066-C|5 

TYR83-C , 




FIG. 8: Radial distribution functions between surface residues 
of PC and oxygens of water. The upper panel shows ionized 
residues and the lower panel refers to non-polar residues. The 
legends in the figure list: aspartic acid (ASP), the probe atom 
is the oxygen at the 1st 5 position; lysine (LYS) , the probe 
atom is nitrogen at the ( position; proline (PRO), the probe 
atom is the /3 carbon; tyrosine (TYR) , with the first e carbon 
as the probe atom. 



the Ox state. Three different charging schemes have been 
used and compared to MD simulations (Tableland Sec. 
HV]) . In the following we will discuss the results relevant 
to charging scheme II only, which are also visualized in 

Fig.m 

The NRFT calculations listed in Table HIl and shown in 
Fig. [7] have been done by using Eqs. ^ and ([9]) in which 
the electric field of PC in Ox state was used for Eo(k). 
The close diamonds in Fig. [7] refer to the total solvent re- 
sponse, while open diamonds represent the polarization 
response only [xp in Eq. ^]. Two interesting observa- 
tions result from examining Fig. [7] (i) a close proximity 



of the NRFT result to the standard (vdW) continuum 
calculation and (ii) a good agreement between the po- 
larization portion of the NRFT calculations and MD re- 
sults. The continuum electrostatics does not reproduce 
the slope of the average energy as we also discuss below 
in relation to the redox entropy. 

In order to understand the origin of the close agree- 
ment between MD and the polarization component of 
the solvent response, one needs to recall what comes to 
the calculation of the polarization and density compo- 
nents of the solvation free energy. The polarization re- 
sponse is calculated by assuming that the only influence 
of the solute on the polarization field is to exclude it 
from the solute volume represented by the step function 
6'o(r) equal to one inside the solute and zero otherwise 
[Eq. mi)]. The density component corrects this result 
by taking into account the inhomogeneous density pro- 
file formed at the surface of a hard-wall solute. In dense 
liquids, such a profile is characterized by a sharp peak of 
the radial distribution function in the first solvation shell 
of the solute. Correspondingly, reflecting the belief that 
the short-range structure of liquids is primarily deter- 
mined by rcpulsionSfSS the density structure factor S{k) 
in Eq. (PT|) was taken in our calculations from the Percus- 
Ycvick solution for hard spheres 1^ The close proximity of 
the full NRFT calculation to the standard DelPhi/vdW 
algorithm (Fig. [7| illustrates the fact that the common 
parametrization of the atomic radii is based on the ex- 
perience learned for hydration of small ions with tightly 
bound first solvation shell. In the present algorithm, this 
physics is accommodated by the density component of 
the solvation free energy. 

The structure of water at the protein surface is quite 
different from what is normally obtained by inserting a 
small solute in a molecular solvent. The structure is het- 
erogeneous including islands of highly structured water 
around polar and ionized residues combined with much 
softer density profile at the hydrophobic patches. This 
reality is illustrated in Fig. [5] which shows pair distri- 
bution functions between ionized and non-polar residues 
and water's oxygens. While the distribution functions 
of ionized residues are reminiscent of the structures typ- 
ically observed around small solutes in dense solvents, 
the water structure around non-polar residues is quite 
different: there is no first-shell peak and water interface 
is shifted by ~ 1 A, in accord with simulations of nano- 
scale hydrophobic solutes 1^ 

The stronger attraction of the surface water molecules 
to the bulk than to a non-polar hydrophobic patch of 
the protein (cavity expulsion potentia l^^'^^ ) results in a 
weak dewetting of the surfacei^ with the density at the 
interface lower than in the bulk (Fig. [8]). Since there are 
only a few charged residues on the protein surface, the 
average surface structure is closer to a step-wise cut-off 
introduced in the polarization component of the response 
function than to a structured liquid at the surface of a 
small polar/ionic solute;^ This observation explains a 
good agreement between MD and polarization calcula- 
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FIG. 9: Solvation Gibbs energy from MD simulations (closed 
circles), NRFT calculations (closed diamonds), continuum 
DelPhi calculation with the standard cavity definition (up- 
triangles), and the cavity surface augmented by the solvent 
radius (t/2 (down triangles). The dashed lines are linear re- 
gressions through the points. 
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FIG. 10: Reorganization energy of PC vs temperature cal- 
culated from MD simulations (closed circles), from NRFT 
(diamonds), and from dielectric continuum using solvent- 
accessible cavity definition (triangles). The dashed lines are 
linear regressions through the points. The filled diamonds 
refer to the full NRFT calculation and the open diamonds 
denote the polarization response only. 



tions of the solvation thermodynamics in this paper as 
well as an equally impressive agreement virith the sim- 
ulations obtained in our previous calculations of charge 
transfer across a polypeptide bridgcj^^ 



C. Solvent Gibbs and reorganization free energies 

The results of calculations of the redox solvation en- 
ergy and solvent reorganization energy [Eqs. (j4ip and 
(|42|l ] using different levels of the theory are listed in Table 
mil Redox and reorganization entropies are given in Ta- 
ble IIVI In addition, the temperature dependence of Ajig 
and As are visualized in Figs. [QlandfTOl The differences 
in the theoretical results arise from the different level of 
structural solvent information incorporated in each for- 
malism. For completeness, we have also listed in Table 
mil the NRFT calculations with the solvent parameters 
of ambient water. 

The continuum electrostatics gain access to the solva- 
tion entropy through the temperature variation of the 
solvent dielectric constant. Therefore, since the dielec- 
tric constant commonly decreases with heating, the sol- 
vent becomes effectively less polar and the solvation free 
energy of a charge distribution increases, i.e. becomes 
less negative. If one considers redox species positively 
charged in both redox states, the oxidized state carries 
a larger charge and hence the difference of Red and Ox 
solvation energies has a positive value decreasing with 
increasing temperature. The redox entropy in Eq. (|40p is 
then positive as is typically observed for simple inorganic 
ionSi ^^i^^ By the same arguments, the species carrying 
negative charge in both redox states should have a neg- 
ative redox entropy, which is the case for the negatively 
charged PC in our calculations. 

Despite the right sign of the redox entropy, the mag- 
nitudes of both the Gibbs solvation energy and the en- 
tropy are markedly different in continuum and micro- 
scopic/simulation approaches: Afis is higher in the stan- 
dard implementation of DelPhi (vdW cavity) than the 
NRFT value by a factor of 1.5 while the redox entropy is 




FIG. 11: Pair distribution function between oxygen of water 
and Cu of PC in the reduced (Red) and oxidized (Ox) sates. 



lower by a factor of ten. The use of the solvent-accessible 
cavity brings the value of A/x^ in a closer proximity to 
the NRFT, but the redox entropy is lowered even more 
(Table [IV|. The magnitude of Afj.^ from NRFT is sig- 
nificantly higher than from MD even if the density com- 
ponent is subtracted from the total response. This ini- 
tially comes a bit of surprise given a good agreement be- 
tween the polarization-NRFT and MD values of (Vbs)ox 
reported in Table Hi] and Fig. [71 We do not currently have 
a good explanation of this disagreement (see Discussion 
below) . 

Concerning the reorganization energy calculations, the 
standard DelPhi/vdW algorithm gives Xg three times 
larger than NRFT and almost an order of magnitude 
larger than the MD simulations (Table IIVP . The origin 
of large Ag in the standard (vdW) implementation of the 
continuum model is in placing highly polar dielectric into 
the small pocket near copper which water molecules do 
not visit in MD simulations. Figure [11] shows the pair 
distribution function between the Cu ion and oxygen of 
water testifying to the fact that water never comes to 
Cu closer than 6 A and the maximum of the first solva- 
tion shell appears at 6.7 A. The addition of the solvent 
radius to the cavity corrects this error reducing the reor- 
ganization entropy to the level of 0.64 eV consistent with 
the polarization component of the NRFT (0.75 eV, Table 
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TABLE IV: Redox solvation free energy Afis and redox entropy Ass in for the Red/Ox states of PC. Also listed are the 
reorganization energy and reorganization entropy, sx — —dX/dT. All energies are in eV and entropies are in meV/K, T = 300 
K. 



Method 



Afis 



ASs 



\s 



S\ 



DelPhi with vdW cavity 

DelPhi with solvent-accessible cavity 

Non-local polarization response functions" 

Molecular Dynamics' 

Experiment 



9.92 


-1.25 


3.62 


0.97 


7.12 


-0.45 


0.64 


0.005 


6.01(-1.33) 


-10.5(4.8) 


1.30(0.55) 


3.2(2.0) 


2.81 


-7.40 
-0.4" 


0.54 


10.7 



-1.4" 



"The numbers in the parentheses indicate contributions to the 
solvation free energy and entropy from density fluctuations. 
''MD results refer to hnear fits through the simulation points. 
=Ref.[23 
'^Ref.m 
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FIG. 12: Reorganization energy of charge-transfer transition 
in p-nitroaniline dissolved in SPC/E water. The results are 
obtained by MD simulations^^ (closed circles) and NRFT cal- 
culations (diamonds) . Closed and open diamonds refer to the 
full and polarization response, respectively. Triangles denote 
half of the Stokes shift from MD simulations; in the LRA, 
Aust/2 = As. The dielectric constants of SPC/E water at 
different temperatures, required for the NRFT input, were 
taken from MD simulations.— 



ences between our calculations and MD simulations. In 
order to illustrate this point, we show in Fig.[T2]the com- 
parison of the NRFT calculations with our recent MD 
results for solvation of a small charge-transfer molecule 
p-nitroaniline in SPC/E water.— The partial charges for 
this molecule, used in the calculations were tabulated 
in Ref. [93. The comparison with the NRFT method is 
somewhat complicated in this case by nonlinear solvation 
effects seen in the deviation of the half of the Stokes shift 
from the reorganization energy. However, the NRFT cal- 
culation with the density component included is in rea- 
sonable agreement with the simulation data for both the 
reorganization energy and entropy, while the polarization 
response alone clearly underestimates the reorganization 
energy. 



VI. DISCUSSION 



The combination of the absolute values of the reorga- 
nization energy with the reorganization entropies clearly 
indicates that re-scaling of the dielectric cavity, often em- 
ployed in various continuum formulations, does not solve 
the solvation problem. Entropies calculated by includ- 
ing the effect of density fluctuations on the solvent re- 
sponse are generally in better agreement with MD sim- 
ulations than the results obtained from polarization re- 
sponse only. All solvation free energies obtained from 
such calculations, however, significantly exceed the sim- 
ulation results. The problem lies in the assumed dense 
structure of the liquid around the solute, which is obvi- 
ously not realized for hydrated protein. Capturing den- 
sity fluctuations is essential for the correct calculations 
of the entropies, but the effective density should be ad- 
justed to that of the structurally loose hydration layer at 
the protein surface. 

One might argue that the failure to reproduce the re- 
sults of the simulations by the NRFT calculations can 
be traced back to the deficiencies of the model and not 
to the specific structure created by the protein in water. 
We believe that this objection cannot explain the differ- 



The formalism developed here is based on the LRA 
suggesting that knowing the variance of electrostatic fluc- 
tuations around a neutral repulsive solute is sufficient to 
calculate the solvation chemical potential [Eq. ^]. One 
can arrive at Eq. ([2]) from the following simple considera- 
tions. The chemical potential of solvation can be written 
as the integral over the magnitude of the solute-solvent 
electrostatic interaction e = Vqs as follows 



-/3mos 



deP{e,(i)e 



-/3e 



(43) 



Here the probability density of reaching the value e is ob- 
tained by taking the statistical average over the reference 
Hamiltonian Hq which includes all the interaction poten- 
tials in the system except the solute-solvent electrostatic 
potential Vqs'- 



P{e, 13) = Qo-i j 5{e- Vq,) e-'^^^dF, 



(44) 



where Qq = Jexp[— /3iJo]rfr and F denotes the phase 
space volume^ Since the solute-solvent component of Hq 
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includes mostly isotropic short-range interactions, one constant background potential $ Eq. (|4ip modifies to 



can put (Vois}o = 0- In addition, the Gaussian approxi- 
mation can be applied to P(e, P) 



P{e, (i) oc exp 



e 
2ai 



with 



{{5Vosf)o- 



(45) 



(46) 



Combining Eqs. (j43p and (|45p . one immediately arrives 
at Eq. (O. 

One of the important lessons of Eqs. (|43| and (144]) 
is that all the thermodynamic information required for 
the solvation problem is contained in the distribution of 
electrostatic potential fluctuations around a "non-polar" 
solute in which all the electrostatic interactions with the 
solvent have been switched off. For large solutes, the 
spectrum of electrostatic fluctuations around the repul- 
sive solute core will ultimately determine the thermody- 
namics of solvation. It is important to stress that Eqs. 
p3)l and ((44|) are exact and this statement is not limited 
to linear solvation only. 

The reasoning outlined above works well for small and 
medium-size solutes, as is exemplified in Fig. 1121 How- 
ever, the application of the same procedure to protein 
solvation studied in this paper has encountered some dif- 
ficulties suggesting that methods developed over several 
decades and successfully applied to solvation of small so- 
lutes in dense polar solvents are probably not directly 
transferable without significant modifications to meso- 
scopic hydration of proteins. The application of our al- 
gorithm to the calculation of the average solute-solvent 
interaction energy turned out to be quite successful when 
the density profile of water around the protein is approxi- 
mated by a step function (Fig. [7|) . The calculations agree 
with MD within simulation uncertainties (< 5%). As 
mentioned above, this comes as a result of averaging be- 
tween tight and loose water structures at the protein sur- 
face. The application of the same procedure to solvation 
of the difference charges of the active site (reorganization 
energy. Fig. [T0|) was less successful, but the agreement is 
probably still acceptable, in particular at lowest temper- 
atures. 

Where the calculations and simulations come in sig- 
nificant disagreement is for the free energy of solvation 
obtained in simulations as the mean of two vertical tran- 
sition energies [Eq. (HI])]. Some recent simulations of cav- 
ities in force-field water modelsiSSiiiSi and of uncharged 
proteiniS^ have suggested a possible origin of this effect. 
It was found that water structured at the protein sur- 
face creates a positive potential within an uncharged cav- 
ity/protein. In terms of Eq. (|45|l this implies a constant 
shift of the solute-solvent energy e -^ e — eo. Since we 
have assumed random orientations of water around an 
uncharged solute, as is indeed the case for small solutes^^ 
our calculations do not include the effect of a positive 
background potential and include only changes of the po- 
tential in response to protein's charges. In the case of a 



A/i, = ((AFo,,)ox + (Ayo.)Rcd) /2 - Ag$, (47) 

where Ag = q^ied — (Zox = —1- AshbaughiSS reported 
an average positive potential of about e$ ~ 9 kcal/mol 
for cavities in SPC water comparable in size with PC. In 
terms of our calculations, it amounts a positive shift of 
the simulation data by about 0.4 eV which will increase 
the current gap of about 3.2 eV (300 K) between the MD 
and NRFT. 

The origin in the difference in solvation free energies 
between calculations and MD simulations might be re- 
lated to the weakly dewetted water density profile near 
the active site. This would imply that some of the prop- 
erties of water structure around large cavities expelled 
by the protein from its volume are quite different from 
the common experience gained with small solutes. This 
qualitative difference between small-size and large-size 
solvation has recently gained appreciation for hydropho- 
bic solvation^ as we discuss next. 

The Gaussian model of hydrophobic solvation goes 
back to the Pratt-Chandler theory of hydrophobicityiS^ 
recently extended by Pratt and co-workers J^^j^S^ The 
formulation of the theory of hydrophobic solvation fol- 
lows a path similar to the one outlined in Eqs. (|43|) - 
(|46ll asking what is the free energy /xn needed to sol- 
vate a solute of volume VLq. It is given by the Gaussian 
probability22ii2i 



/3mo 



2xo 



(48) 



with the fluctuation of the number of solvent particles in 
volume ilo 

Xn = {{5Nf)n = P^o + p" I drdr'hU\r - r'|). (49) 

Jna 

In Eq. (|49|) . hss{r) is the pair correlation function of the 
homogeneous solvent. 

The Gaussian probability of electron transfer carries a 
close similarity with Eq. (|48l) giving the activation free 
energy as 



Pp' 



y2 
act „ ^0 

~ 2ar 



(50) 



where Xq is the average vertical donor-acceptor energy 
gap and a^ is the variance of the solute-solvent interac- 
tion potential when Vqs — AVqs is used in Eqs. (|44|) and 
(|45|) . Not surprisingly, the structure of the equation for 
a^ resembles Eq. ((49|) : 



<y', = P J VoMl)dr, 

+ p^ J Vo,il)Vos{2)hUl, 2)dridr2 (51) 
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The last sumniand in this equation represents the density 
component of the response since it transforms into the 
k-space integral with the density structure factor of the 
liquid [Eq. dH])]. 

Equation (j5ip carries a close resemblance with the 
Pratt-Chandlcr theory of hydrophobicityJ^ The sum 
of the first and the third terms is the average of the 
squared solute-solvent potential V^^s(l) over the solute- 
solvent density profileiS^ 

p(r)=p(l-fco,s(r))+p2 /co.(r')/is.(|r-r'|)dr' (52) 



in which the solute-solvent direct correlation function 
cos(r) is replaced by its lowest density expansion^S 
cos(r) ~ — ^o(r) (^o(r) is one inside the solute and 
zero otherwise). The observation that the density pro- 
file around solutes of size > 1 nm can be approximated 
by a step function (Fig. [TI] also see Fig. 3 in Ref. l97l ) 
amounts to neglecting the second term in Eq. (j52p and, 
correspondingly, the density components Xd in the sol- 
vent response function [Eq. p^ ]. 

The fact that the spectrum of electrostatic poten- 
tial fluctuations around a non-polar solute gives com- 
plete information about polar solvation brings this latter 
problem in close relation to the problem of hydropho- 
bic solvation. It was realized in recent years that hy- 
drophobic solvation of small and large solutes are quali- 
tatively different.— i2iii2£ Solvation character changes at 
the critical size of ~ 1 nm from entropy-dominated 
solvation of small solutes (Gaussian statisticsiSl) to 
enthalpy-dominated solvation of large solutes driven by 
the creation of the solute-solvent interface (non-Gaussian 
statisticsiSl). The solvent interface around large so- 
lutes involves partial dewetting, strongly sensitive to the 
strength of solute-solvent attractions^^ and the appear- 
ance of large-size interfacial density fluctuations i^ In 
case of protein solvation, surface water creates a nonzero 
potential around uncharged proteins discussed above, 
while density fluctuations lead to complex protein-solvent 
dynamicsii and arc probably connected to "slaving" of 
the protein dynamics by the solvent ji^ 

These new features observed for hydrophobic solvation 
at the nano-scale will affect the thermodynamics of po- 
lar solvation. It is currently not clear how the cross- 
over from the Gaussian regime of small solutes to the 
interface-dominated regime of large solutes will trans- 
late to the problem of electrostatic solvation. The cal- 
culations and simulations performed in this study gave 
some insights into the kind of problems which the the- 
ory needs to address. It is clear that the inhomogc- 
neous nature of the protein surface requires algorithms 
for the local density profile^ or the solvent-accessible 
surfacoiSS to be a part of a qualitative solvation theory. 



The current formulation can be improved by using a lo- 
cal density approximation as, for instance, applied by 
Ramirez and Borgis.— In this approach, the bulk dipo- 
lar density y [Eq. ([6])] is replaced by the local dipolar 
density y{r) = (47r/9)/3TO'2p(r) -I- (47r/3)ap(r) with the 
local density profile calculated from, for instance, the 
Lum- Chandler- Weeks theoryiiS^ Our computational al- 
gorithm will then require the Fourier transform of the 
field Eo(r)p(r)/p instead of Eo(r)(l — 0o(r)) in the cur- 
rent implementation. In addition, density fluctuations 
of the interfacial region need to be addressed since the 
mean position of the interface has no physical significance 
in the presence of large-amplitude fluctuationsjsii What 
is also clear is that the density fluctuations of the interfa- 
cial region present a nuclear mode, largely diminished for 
solvation of small solutes, which grows in significance for 
solvation of biopolymers. Future theoretical development 
needs to address this physical reality. 

Both NRFT and MD simulations predict a substan- 
tial temperature dependence of the redox potential (ca. 
~ 5 — 7 mV/K). This magnitude of the temperature vari- 
ation is prohibitively high since a temperature change of 
~ 15 K would shift the redox potential by about 100 
mV potentially terminating many enzymetic redox reac- 
tions. The theoretically predicted redox entropies refer 
to zero ionic strength and the solvent contribution to the 
redox potential. The contribution from the protein to 
the overall redox potential is available from our simula- 
tions, but it depends weakly on temperature contribut- 
ing only ~ —0.9 mV/K to the redox entropy. The com- 
bined solvent/protein entropy is then significantly higher 
than redox entropies experimentally observed in buffered 
solutionsSii^ (Table ITV]) . This difference raises the ques- 
tion of the role of the ionic atmosphere in stabilizing the 
redox potential of mettalloproteins. The existing experi- 
mental evidence for small redox molecule a^^'^i^^^ and sim- 
ulations of metalloproteinsS^ all indicate a substantial ef- 
fect of the ionic atmosphere on the redox entropy, which 
might compensate for the large entropy due to water and 
protein. Since simulations have little to say about the 
ionic strength effects, laboratory measurements of redox 
entropy at different buffer concentrations are required to 
shed more light on this problem. 
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